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ABSTRACT 

We derive and present a fast and accurate solution of the initial value problem for 
Keplerian motion in universal variables that does not use the Stumpff series. We find 
that it performs better than methods based on the Stumpff series. 
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1 INTRODUCTION 


IWisdom fc HolmanI (|l99lh introduced a symplectic mapping 
method for the rapid simulation of the n-planet problem (n 
planets plus a massive central body). The method splits the 
Hamiltonian for the n-planet problem into Kepler Hamilto¬ 
nians and an interaction Hamiltonian, each of which may be 
efficiently solved. The evolution of the n-planet problem is 
obtained by interleaving the elementary pieces. The rapid 
and accurate solution of the Kepler initial value problem 
is an essential part of the method. The Wisdom-Holman 
method and its variations has been widely adopted for solar 

system dynamics i nvestigations. _ 

The method of lWisdo m H olman Jl 99 ih evolved from 
the mapping method of IWisdoml (|l982l ): it relies on the 
averaging principle to introduce Dirac delta functions into 
the Hamiltonian. An alternate approach, which leads to the 
same algorithm, is that of symplectic integration, which 
uses the algebra of Lie series to approximate the local evo¬ 
lution to some order in the stepsize by interleaving the 
evolution of the pieces. An advantage of the mapping ap¬ 
proach is that the stability of the method can be analyzed in 
terms of the overlap of “stepsize resona nces,” which can be 
read off the delta function Hamiltonian ([Wisdom HolmanI 
I 1992 I ). Another advantage of the mapping approach is that 
perturbation theory can be used to improve the method 
by eliminating the high-frequency terms introduced by 
the delta functions, leading to th e “symplectic corrector” 
([wisdom. Holman Touma Il996h . An advantage of the 
symplectic integration approach is its a lgebraic simplicity. 

_ Jacobi coordinates were used in I Wisdom HolmanI 

([l99lh to eliminate the center of mass freedom and to effect 
the separati on of the Hamiltonian into Kep lerian and inter¬ 
action parts, [^uma fc WisdomI ([l993l . [l993 ) used a different 
splitti ng, making use of the ca nonical heliocentric coordi¬ 
nates. Ibevison DuncanI ([l994h used the Wisdom-Holman 
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method in Jacobi coordinates, patched together with a non- 
symplectic method for handling close encounters among the 
bodies. They distribu ted their program, call e d RMV S3, in 
their SWIFT package. IPuncan. Levison &: Leel ([ 19981 ) used 
the canonical heliocentric variables with a slightly differ¬ 
ent splitting; they call their method the “democratic he¬ 
liocentric” method. They handle close encounters by recur¬ 
sively subdividing the step in a s ymplectic manner . Their 
program is distributed as SYMBA. IChambersI ([l999l ) intro¬ 
duced another method based on the democratic heliocen¬ 
tric splitting, with a symplectic transition to an ordinary 
numerical integration method (the Bulirsch-Stoer method) 
for close encounters . His p rogram is distributed as MERCURY. 
Ibevison Duncai] ([2000 ) noticed that methods based on 
the democratic heliocentric splitting are unstable for large 
eccentricities, and modihed their SYMBA program to numeri¬ 
cally integrate “close encounters with the Sun.” Their mod¬ 
ified program is called “modified SYMBA.” An essential el¬ 
ement of all of these variations on the Wisdom-Holman 
method is that one must solve the Kepler initial value prob¬ 
lem: given the position and velocity at one time, the task 
is to find the position and velocity at a different time (dis¬ 
placed by a timestep that can be either positive or negative). 

A recent contribution using the Lie series approach is 
iHernandez Bertschinge3 ([20151 ). They developed a sym¬ 
plectic integrator for the collisional n-body problem. This 
method also relies on the rapid and accurate solution of the 
Kepler initial value problem. 

We have found that the solution of the initial value 
problem may be efficiently and accurately carried out in 
universal variables. The resulting formulation and program 
work for all orbits, whether they are elliptic, parabolic, or 
hyperbolic. The t r aditio nal presentation of univeral vari¬ 
ables (e.g. Id an ^ (I 1992 I 1 1 makes use of Stumpff functions 
and calculates them using their series representation. Ar¬ 
gument four-folding is used to bring the function argument 
into an interval near zero so that the series converge rapidly 
enough. However, all of the functions can be represented in 
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2 J. Wisdom and D.M. Hernandez 


terms of ordinary trigonometric functions. But, a straight¬ 
forward computation using this representation has unpleas¬ 
antly large error, so the series and argument four-folding 
appear to be necessary. We have been able to circumvent 
this problem by reexpressing the functions in a numerically 
well defined way. The builtin trigonometric functions are 
fast and accurate, and our solution of the initial value prob¬ 
lem compares favorably with solutions based on the Stumpff 
series. 


2 KEPLER PROBLEM 


The Hamiltonian for the Kepler problem is: 


2m r 


(1) 


where p is the magnitude of p and r is the magnitude 
of X. The constants fi and m depend on the context in 
which the Kepler problem is found. For example, if x = 
X 2 — xi is the relative coordinate in the two-body problem, 
then we would be led to take m to be the reduced mass 
(1/mi + l/m 2 )~^ and ji — Gmim 2 - In the n-body problem, 
the const ants may be different, dep ending on the formula¬ 
tion fe.g. IWisdom HolmanI (|l99lh b 

Hamilton’s equations may be written as a second order 
system: 


X = 


/cx 
^ ’ 


( 2 ) 


where the “Kepler constant” k — p/m. The evolution de¬ 
pends on the constants only through this combination. We 
use the dot notation for derivative with respect to time. 


3 DERIVATION OF THE SOLUTION 


A basic ref erence for the solution of the Kepler initial value 
problem is IPanbvl (|l992h . We refer the reader to that pre¬ 
sentation of universal variables and the Stumpff series for 
background on the series approach. We follow some aspects 
of that presentation here. Our derivation is self-contained. 

We introduce a new independent variable s satisfying 


ds 1 
dt r 

The equation of motion, Eq. becomes 


(3) 


x" - Lx' + -X = 0, (4) 

j, j, 

where prime indicates differentiation with respect to s. 

The Hamiltonian is conserved, since there is no explicit 
time dependence. It is conventional to write the conserved 
quantity as 



a 


2k . . 

-X X. 

r 


(5) 


The constant /3 is positive for elliptic motion, for which a 
is the semimajor axis, (3 is negative for hyperbolic motion, 
and /3 is zero for parabolic motion. 

Expressing /3 in terms of derivatives with respect to s 
yields: 


/? = 



( 6 ) 


Differentiating this expression with respect to s and using 
the derivative of the equation of motion, Eq. yields: 

:k"' + /3x' = 0. (7) 

Next we introduce the / and g functions: 

X = /xo+^Vo, 

V = /xo + ^vo, (8) 

where xo is the initial position, and vo is the initial time 
rate of change of position (the initial velocity). Substituting 
this into Eq. ©, we find 

f' + Pf = 0 

g'” + Pg' = 0, (9) 

using the independence of initial position and velocity. The 
equations are satisfied by an offset simple harmonic oscilla¬ 
tion in s. 

To determine the initial values for the derivatives, we 
start with Eq.@, to hnd 

/ + ^/ = 0 

3+—S = 0. (10) 

Erom Eq. m, the initial values satisfy /o = = 1 and 

/o = Po = 0. Then 

/o = rofo = 0 
9o = Togo = ro 

fo = ^o/o + rorofo = - — 

ro 

9o = rlfo + rorofo = roro- (11) 


It is convenient to dehne solutions of Eq. Q in terms 
of functions (s) satisfying 


Gf(s) = folds)- 

(12) 

We start the ladder with 


Go(s) = cos(y^s), 

(13) 

for /3 > 0, and 


Go(s) = cosh(y^s), 

(14) 

for /3 < 0. We can take 


sin(^/ps) 

Ods)- ^ , 

(15) 

for /3 > 0, and 


sinh(V^s) 

’ 

(16) 

for /3 < 0. We can then take 


G^s) = (1 - cos{^/ps))/P, 

(17) 

for /3 > 0, and 


Ols) = (1 - cosh(y^s))//3, 

(18) 

for /3 < 0. Onward, 

G^s) = {s-sm{^/ps)/-/p)/p 


= {s-Gls))/p. 

(19) 
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Universal Kepler Solver 3 


for /3 > 0, and 

G§(s) = {s - sinh{^/^)/^/^)/P 

= (s-Gf(s))//3, (20) 

for /3 < 0. We could go on, but this is all we need here. 
We have chosen the constants so that (0) = 0 for z > 0. 
Notice that for z < 3 these functions satisfy 

(Gf(s))'" + /3(Gf(s))' = 0. (21) 

Now we can use these functions to find solutions for / 
and g. Let 

f{s) — AfGl{s) + BfG2{s) + Cf. (22) 

The condition that /(O) = 1 implies Gf — 1. Then we form 
the derivative 

f\s) = AfG^Us) + B}Gl{s). (23) 

The condition that f'{G) — 0 implies Af = 0. The next 
derivative is 

fis) = BfG^Us), (24) 

using the fact that Af = 0. The condition that f"{0) = 
—k/ro implies that Bf = —k/ro. Putting it together we find 

/(s) = 1 - (fc/ro)G§(s). (25) 

Similarly, we let 

g(s) = AgGi{s)BgG^is)Gg. (26) 

We find 

g{s) = roG^{s) + roroG^{s). (27) 

From these we find / and g by using the relation dsidt — 
1/r to hnd 

f{s) = -(fc/(rro))Gf(s), 

g{s) = (ro/r)(G^(s)+roGf(s)). (28) 

But we still have to find s\ 

The relation between t and s is 



h — t — to — / r{s)ds. 

Jo 

(29) 

Let’s express r(s) in terms of Cf(s). First, 



r'o = roro = xo • xq. 

(30) 

After a 

small reduction, we hnd 



0 

II 

1 

0 

(31) 

and 

C + pr'o = 0 . 

(32) 

So r(s) 

may be expressed as 



r(s) = ro + TqGi (s) + VqG^ (s). 

(33) 

Using Eqs. (l30l) and (l3T1) 



r(s) = roCo(s) + roroCf (s) + kG^is). 

(34) 

The properties of the Cf allow an immediate integration to 
find 


h — roC^(s) + TQT 0 G 2 + kG^i^s^. 

(35) 


This implicit equation for s is our “Kepler equation.” We 


can solve this by any of the standard methods: Newton, Hal¬ 
ley, Laguerre, bisection, and so on. In practice, we first try 
Newton’s method. If this fails, we try the Laguerre-Conway 
method. If this fails, we recursively subdivide the step. 

For the hyperbolic case our solution is as follows. For 
small stepsizes the equation is well approximated by a 
cubic equation. We take the real solution of this cubic 
equation as our initial guess in Newton’s method. If New¬ 
ton’s method does not converge, then we use the Laguerre- 
Conway method. If this fails we recursively subdivide the 
step. 

Eq. (IMl) can be used to rewrite the equation for g. We 

find 

p(s) = 1 - (fc/r)G§(s). (36) 

The expressions for the parabolic case can be found by 
taking a limit as /3 goes to zero. As it turns out the Kepler 
equation becomes a cubic equation in s and so can be solved 
without iteration. 


4 NUMERICAL REFINEMENT 

In order to have expressions that are well defined numeri¬ 
cally we have to do a little more work. 

Examining Eqs. (HZl and m, we see that for small ar¬ 
guments of the cosine and hyperbolic cosine functions there 
will be cancellation and loss of precision. But this problem 
can be fixed by using the half-angle formulas. 

Assuming /3 > 0, let 

S 2 = sin(v^s/2) 

C 2 = cos(v^s/2), (37) 

then we can write 

Cf(s) = 2s2C2/y/;d 

G^{s) = 2 S 2 S 2//3 

Gsis) = (s-Cf(s))//3 

Cg(s) = l-/3C^(s). (38) 

The only expression that might be of concern is the expres¬ 
sion for C 3 (s), but, in practice, it seems to not be a problem. 
Similar expressions can be derived for the /3 < 0 case. Such 
changes are made throughout the program. 


5 NUMERICAL EXPLORATION 


We compare here our method and program (universal. c), 
to two universal variable Kepler stepper programs that use 
Stumpff series. The program drift.one.f was written by 
Harold Levis on and Martin Duncan. It is based on the pre¬ 
sentation in Id an ^ (Il992(l . The program drift.one.f is 
available as part of the SWIFT package. This program is used 
in a nu mber of programs derived f rom the Wisdom-Holman 
method (IWisdom fc Holm^ 199lh . These programs include 
RMVS3 (|Levison &: DuncanI 1 19941 ). also distributed in the 
SWIFT package. The same program i s used to advance the Ke¬ 
pler p roblem in the pr ogram SYMBA ([Puncan. Levison fc Led 
119981 ) , and in MERCURY (|Chamberslll999h . In order to compare 
to our program, written in the C programming language, 
we have translated, line by line, the program drift.one.f 
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to a C version drift.one. c. We have not compared our 
code directly to the fortran program drift_one.f. We also 
compare our program to the Kepler stepper in WHFast 
(|Rein Tamavoll2015l h 

Our test is designed to work for both elliptic and hyper¬ 
bolic orbits. Let T = 27r/n where n = y/k/\a\^. For elliptic 
orbits, T is the orbital period, and n is the mean motion. 
We choose G = (0.0172)^, which approximates the gravita¬ 
tional constant in units of AU, day, and solar mass. We take 
the mass factors to be one; thus the Kepler constant fi/m 
is numerically just G. We use the same Kepler constant in 
the drift.one.c calculations. We take the semimajor axis 
for the elliptic case to be 0.4AU. The eccentricity and step- 
size are varied. The pericentric distance is ^ = a(l — e), and 
the velocity is determined from a through the value of the 
energy. We evolve the Kepler orbit back and forth through 
pericenter, adjusting the phase so that the pericenter is en¬ 
countered with a wide range of phases. The chosen stepsize 
is h; we make use of the auxiliary stepsize h' = 7/1, where 
7 = (\/5 - l)/2 0.618, is the irrational golden mean. We 

start at pericenter with t = 0, and evolve the orbit with 
stepsize h until t > Tl2, i.e. until we have passed apocenter 
(for elliptic orbits). Then we adjust the phase with a pos¬ 
itive step of h' . At this point we start collecting statistics. 
We reverse the timestep and evolve the orbit with timestep 
—h until t < —T12. Then we adjust the phase with a posi¬ 
tive step of h!. We evolve with stepsize h until t > T/2; we 
adjust the phase with a postive step of h\ and so on. We re¬ 
peat this process for 100 pericenter passages. At the end we 
collect statistics again, and compare to the initial statistics. 
This procedure tests all parts of the evolution for a large 
variety of phases. It also tests both positive and negative 
timesteps. 

The results for elliptic orbits for three methods are sum¬ 
marized in Figs. [3l[5l The top plot in each set shows the en¬ 
ergy error as a function of eccentricity e and stepsize h. The 
middle plot shows the sign of the energy error at the end of 
each evolution. We would like this plot to show an intimate 
mixture of positive and negative results, so that an evolu¬ 
tion computed with the Kepler solver does not show any 
tendency to expand or contract . This “bias” plot was intro¬ 
duced by I Rein Tama^ (|2015l b The bottom plot shows an 
estimate of the computing time per Kepler step. The tests 
were run on a 4 GHz iMac, with an Intel i7 chip. All of the 
tests were compiled with C compiler optimizer option -03. 
Comparing Figs. [3]and[31 we see that our code is more ac¬ 
curate, faster, and shows less bias than drift_one.c. Note 
the complicated structure in the bias plot for drift.one.c, 
and the intimate mixture of colors in the bias plots for 
universal.c and the Kepler stepper in WHFast. 

We have computed the average ratio of the time per 
step for the interval 0.001 < h/T < 0.1, which is the range of 
most interest for practical calculations, and find the average 
ratio for drift.one.c relative to universal.c is about 1.9. 
The code universal. c is about twice as fast as drif t_one. c. 
We have also computed the average ratio of the time per 
step for the interval 0.001 < h/T < 0.1 for the Kepler solver 
in WHFast and find that the time per step is about a fac¬ 
tor of 0.79 smaller than our code—it is about 20% faster. 
In Fig. [H histograms of the relative energy error in the el¬ 
liptic case for the three methods are shown. We see that 
drift.one.c is less accurate than both universal.c and 
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Figure 1. These plots histogram the relative error in the elliptic 
case for three methods. The bottom plot is for universal.c, the 
middle plot is for drift_one.c, and the top plot is for the Kepler 
solver in WHFast. 


the Kepler solver in WHFast. Quantitatively, the averages of 
the common logarithm of the absolute value of the relative 
error are: -11.92, -11.74, and -11.09, for universal.c, the 
Kepler solver in WHFast, and drift.one.c, respectively. The 
error of universal. c is, on average, about 50% smaller than 
the Kepler solver in WHFast, and about a factor of 6.5 smaller 
than the error in drift_one.c. 

For the hyperbolic case, we let a = —0.4AU. We start 
at pericenter. The pericentric distance is g = a(l — e); e is 
larger than 1. The initial velocity is determined from the 
energy. The back and forth procedure is the same as in the 
elliptic case. 

The results for hyperbolic orbits for our code, 
un iversal. c, and for dri ft_one. c are shown in Figs. [6] and 
[71 [Rein &: Tamavol (|2015l j did not extensively test the hy¬ 
perbolic case, and the Kepler solver in WHFast performs 
poorly for unbound orbits. Note again that the bias plot 
for our method universal. c shows an intimate mixture, 
whereas the bias plot for drift.one.c shows evidence of 
bias. We find that the average ratio of the time per step of 
drif t_one. c to the time per step of universal. c for the in¬ 
terval 0.001 < h/T < 0.1 is about 1.6; drift.one.c is about 
60% slower than universal.c. 

In Fig. [21 histograms of the relative energy error in 
the hyperbolic case for drift.one.c and universal.c are 
shown. We see that drift.one.c is typically less accurate 
than universal.c. Quantitatively, the averages of the com- 
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Figure 2. These plots histogram the relative error in the hyper¬ 
bolic case for two methods. The bottom plot is for universal.c, 
and the top plot is for drift_one.c. 

mon logarithm of the absolute value of the relative error 
are: -11.72 and -11.03, for universal.c and drift.one.c, 
respectively. The error of universal. c is, on average, about 
a factor of 5 smaller than the error in drift_one.c. 


6 SUMMARY 

We have developed and tested a universal variable solver 
for the Kepler initial value problem. The method eschews 
the use of Stumpff series in favor of ordinary trigonometric 
functions. We have been careful to make sure that all expres¬ 
sions are numerically well defined. We find that our program 
performs better, in terms of accuracy, bias, and speed, than 
a C version of a widely used Kepler solver, which uses the 
Stumpff series. Our code is freely available; contact the au¬ 
thors. 


ACKNOWLEDGEMENTS 

We thank Hanno Rein, Daniel Tamayo, and Edmund 
Bertschinger for helpful discussions. DMH acknowledges 
support by a NSF Graduate Research Fellowship under 
Grant No. 1122374. 


REFERENCES 

Chambers J. E., 1999, MNRAS, 304, 793 
Danby J. M. A., 1992, Fundamentals of celestial mechanics, 
2nd edn. Willmann-Bell, Richmond, Va., U.S.A. 

Duncan M. J., Levison H. F., Lee M. H., 1998, AJ, 116, 
2067 

Hernandez D. M., Bertschinger E., 2015, MNRAS, 452, 
1934 

Levison H. F., Duncan M. J., 1994, Icarus, 108, 18 
Levison H. F., Duncan M. J., 2000, AJ, 120, 2117 
Rein H., Tamayo D., 2015, MNRAS, 452, 376 


Touma J., Wisdom J., 1993, Science, 259, 1294 
Touma J., Wisdom J., 1994, AJ, 107, 1189 
Wisdom J., 1982, AJ, 87, 577 
Wisdom J., Holman M., 1991, AJ, 102, 1528 
Wisdom J., Holman M., 1992, AJ, 104, 2022 
Wisdom J., Holman M., Touma J., 1996, Fields Institute 
Communications, Vol. 10, p. 217, 10, 217 

This paper has been typeset from a T^lX/ DT^]X file prepared 
by the author. 


© 2015 RAS, MNRAS 000, [HS] 












6 J. Wisdom and D.M. Hernandez 



log,oih/T) 

Figure 3. The summary diagrams in the elliptic case for universal.c are plotted. The top plot shows the relative energy error as a 
function of eccentricity e and stepsize h. The color bar shows the common logarithm of the magnitude. The middle plot illustrates bias 
in the calculation, plotting 1 if the energy error is positive, -1 if the energy error is negative, and 0 for zero energy error. The bottom 
plot shows the computing time in microseconds per call to the Kepler stepper. 
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Figure 4. The summary diagrams in the elliptic case for drif t_one. c are plotted. The details are the same as in Fig. |3] 
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Figure 5. The summary diagrams in the elliptic case for the Kepler solver in WHFast are plotted. The details are the same as i: 


Fig. El 
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Figure 6. The summary diagrams in the hyperbolic case for universal.c are plotted. The details are the same as in Fig. |3] 
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Figure 7. The summary diagrams in the hyperbolic case for drift_one. c are plotted. The details are the same as in Fig. |3l 
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